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PERI-TYRRHENIAN PHYLOGEOGRAPHY IN THE LAND SNAIL 
SOLATOPUPA GUIDONI (PULMONATA) 
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ABSTRACT 


The land snail Solatopupa guidoni has an interesting disjunct distribution in the peri-Tyrrhe- 
nian area, with a few scattered populations in Sardinia, Corsica, and Elba Island. The species 
is strictly bound to limestone and has poor dispersal capacity. Here, we used sequences of 
two mitochondrial genes, cytochrome oxidase | (COI) and 16S rRNA (16S), to reconstruct 
the species phylogeography. The study is based on nine of the 13 populations reported in 
the literature. Phylogenetic methods and nested clade analysis failed to retrieve reciprocally 
monophyletic lineages on each island. Sardinian and Elba Island populations were found to 
be embedded within those of Corsica. These results contrast with previous phylogeographic 
studies of other organisms with similar distribution patterns, showing major phylogeographic 
breaks between islands. Molecular clock estimates suggest that species diversification oc- 
curred less than 2.5 Myr ago, in line with a previous study that dated the split of S. guidoni 
from its continental sister species at 3-5 Myr. Corsican populations can be tentatively sorted 
into two haplo-groups. A first group includes haplotypes found in central inland Corsica only, 
whereas haplotypes found in northern coastal locations cluster in a second haplo-group. The 
latter also includes Sardinian and Elba Island populations. Molecular data indicate that a 
combination of vicariance and dispersal events provide the best explanation of the distribu- 
tion of this species. 
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INTRODUCTION 


The genus Solatopupa Pilsbry, 1917, com- 
prises a small group of six land snail species 
with a peri-Tyrrhenian distribution, ranging 
from northern Spain to central Italy, including 
Sardinia, Corsica, and Elba Island. All but one 
species (S. cianensis) ascribed to the genus 
are strictly bound to limestone. These charac- 
teristics, coupled with the general low potential 
for long-distance dispersal of land snails, make 
Solatopupa a reliable model for biogeographic 
studies. The distribution of the genus has been 
traditionally explained by vicariance due to 
paleogeographic events that caused the de- 
tachment of the Sardinia-Corsica microplate 
from the Iberian plate in the early Oligocene, 
its subsequent rotation anti-clockwise towards 
the Italian peninsula (early to middle Miocene), 
and finally its separation into two islands (Alva- 
rez, 1972: Bellon et al., 1977; Boccaletti et al., 


1990; Carmignani et al., 1995; Gelabert et al., 
2002; Bosellini, 2005). The interaction of the 
microplate with the emerging Apennine chain 
during the Pliocene has been invoked to explain 
the occurrence of a single species, Solatopupa 
guidoni (Caziot, 1902), in Sardinia, Corsica, 
and Elba Island (Giusti, 1977). 

Ketmaier et al. (2006) produced a phylo- 
genetic hypothesis for the genus Solatopupa 
based on nucleotide sequences of three mito- 
chondrial and one nuclear gene, showing that 
phylogenetic relationships were consistent 
with the paleogeography of the peri- Tyrrhenian 
area. However, contrary to previous reports on 
other taxa with similar peri- Tyrrhenian distribu- 
tions, molecular clock estimates suggested that 
the diversification of So/atopupa coincided with 
the very last phases of detachment of the mi- 
croplate (late Miocene-middle Pliocene) rather 
than with the initial ones (late Oligocene-early 
Miocene). Molecular data also suggested that 
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the relatively widespread S. similis (Bruguière, 
1792), occurring from northern Spain to eastern 
Liguria, Italy, is the continental sister species 
to the exclusively insular S. guidoni. Ketmaier 
et al. (2006) analyzed four populations of S. 
guidoni (two from Corsica, one from Sardinia 
and one from Elba Island) and found a close 
affinity between the Sardinian and Elba Island 
populations, whereas the Corsican populations 
formed a single clade. Previous allozyme data 
on a slightly different set of populations sug- 
gested two genetically distinct lineages in the 
species (Fst = 0.523 over 28 loci, P < 0.05; 
Boato, 1988; at that time, the name S. simonet- 
tae was still used to identify the species, and 
only later was it synonymized with S. guidoni). 
One group comprised two populations from 
inland Corsica and one included two popula- 
tions from northern coastal Corsica and one 
population from Sardinia. 

Solatopupa guidoni has a very scattered 
distribution, with nine populations known from 
Corsica, two from Sardinia and two from Elba 
Island (Fig. 1; Giusti, 1977). Populations are 
strictly associated with calcareous rocks; these 
are of different ages in the case of Corsica and 
Sardinia (Fig. 1). For this study, we consider- 
ably expanded the sampling in terms of popu- 
lations and individuals, and we sampled nine 
out of the 13 known populations and surveyed 
a total of 60 individuals. We could not find any 
individual in four other localities from which the 
species has been reported in the past (Fig. 1). 
We also surveyed other ecologically suitable 
areas (i.e. with limestone), particularly in cen- 
tral and southern Corsica but without success. 
It therefore seems likely that we covered the 
whole range of the species. 

We used fragments of two mitochondrial DNA 
(mtDNA) genes — cytochrome oxidase | (COI) 
and 16S rRNA (16S) — to produce a phylogeo- 
graphic hypothesis for S. guidoni. COI and 16S 
are the most variable genes sequenced by 
Ketmaier et al. (2006); we therefore expected 
to find sufficient genetic variation to clarify 
relationships at population level. We aimed 
to determine (1) whether mtDNA lineages 
occurring on the three islands are recipro- 
cally monophyletic as expected in light of their 
allopatric distribution and the poor dispersal 
capacity of the species. We also want to (2) 
obtain an approximate temporal framework for 
the intraspecific diversification of the species. 
Because the split between S. guidoni and its 
closest continental relative S. similis has been 
dated at 3-5 Myr (Ketmaier et al., 2006), di- 
vergence in S. guidoni should be more recent. 


Finally, we aim to exploit our wider sampling to 
(3) test the hypothesis of Boato (1988) that in 
Corsica S. guidoni would comprise two geneti- 
cally distinct lineages. We discuss our findings 
in light of the paleogeography of the area and 
phylogeographic studies conducted on a vari- 
ety of groups with peri-Tyrrhenian distribution 
(i.e. newts, subterranean terrestrial and aquatic 
isopods, stoneflies, and flightless beetles) 
(Cobolli Sbordoni et al., 1995; Caccone et al., 
1997; Ketmaier et al., 1999, 2000, 2003a, b; 
Fochetti et al., 2004). 


MATERIALS AND METHODS 


Sampling, DNA extraction, PCR Amplification 
and Sequencing 


Five Corsican populations of S. guidoni were 
sampled, and this material was added to the 
four populations (one from Sardinia, two from 
Corsica and one from Elba Island) analyzed 
by Ketmaier et al. (2006). Details on sampling 
localities and sample sizes are given in Table 
1 and Figure 1. Total genomic DNA was ex- 
tracted following Ketmaier et al. (2006). We 
used the same primers and PCR conditions 
as in Ketmaier et al. (2006) to amplify about 
700 and 500 base pairs (bp) of mitochondrial 
DNA (mtDNA) of Cytochrome Oxidase | (COI) 
and 16S rRNA (16S) gene fragments. PCR 
products were sequenced with an automated 
sequencer (Applied Biosystems 373A) follow- 
ing the manufacturer’s protocols. To promote 
accuracy, both strands were sequenced using 
the primers employed for PCR. The sequences 
have been deposited in GenBank (Accession 
No. GU046364—GU046389; Table 1). 


Population Genetics and Phylogenetic Analyses 


Sequences were edited and aligned using 
Sequencher 4.6 (Gene Code Corporation); 
alignments were also checked by eye. The 
suitability of pooling sequence data from the 
two sequenced genes was assessed by the 
incongruence length difference test (Farris et al., 
1994), as implemented in PAUP* 4.0810 (Swof- 
ford, 2002). The test showed that the two regions 
are not phylogenetically incongruent (P = 0.787). 
We further check combinability by using the 
Bayesian approach detailed in Kauff & Lutzoni 
(2002), with a threshold value of posterior prob- 
abilities = 99% for conflicting relationships. None 
of the few conflicting relationships revealed by 
the separate analyses received reciprocal pos- 
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FIG. 1. Schematic maps showing the sampling localities of the Solatopupa 
guidoni populations included in the study and the geology of the region. (a) 
Geographic position of the three islands - Sardinia (S), Corsica (C), and Elba 
Island (E) - where S. guidoni occurs. In (b) and (c) are respectively shown 
sampling localities in Corsica/Sardinia and Elba Island. Three letter codes 
identify different populations and they correspond to those in Table 1. Black dots 
identify localities where the species had been reported but not found during our 
samplings. Downward diagonal shading indicates the distribution of limestone 
in three islands. Different patterned rectangles close to each population refer 
to the different ages of calcareous rocks inhabited by S. guidoni populations 
(this information is not available for Elba Island). Positioning of three islands 
roughly approximates real geographic relationships but sizes of the islands 
are not in proportion to one another (see scale bars). 
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TABLE 1. Locality, latitudinal and longitudinal data (Lat/Log) used in the nested clade analysis, sample 
size (N), and abbreviations for the Solatopupa guidoni populations included in the study. The last two 
columns report GenBank accession numbers. @Allozyme data from Boato (1988); Data from Ketmaier 


et al. (2006). 
GenBank accession numbers 
Locality Lat/Long N Code COl 16S 
Capo Caccia, Sardinia, Italy 40°33’N, 8°9’E 3 GEeEa DQ305079>» DQ305064> 
Corte, Corsica, France 42°18N, O°9E 16 COR GU046364 GU046381 
GU046365 GU046386 
Sabara, Corsica, France 42°25'N, 9°8'E 4 SBRa DQ305080Þ DQ305062> 
Francardo, Corsica, France 42°23 91E 3 FRAa DQ305077> DQ305063> 
Caporalino, Corsica, France 42°22'N,9°11E 5 CAPa GU046366 GU046379 
GU046367 GU046380 
GU046368 GU046385 
GU046369 
Oletta, Corsica, France 42°37N,9°20E 9 OLE GU046370 GU046384 
GU046371 GU046389 
St. Florent, Corsica, France 42°41'N, 9°19E 10 FLOa GU046372 GU046382 
GU046373 GU046383 
GU046374 GU046387 
Sisco, Corsica, France 42°48'N,9°29E 8 SIS GU046375 GU046388 
GU046376 
GU046377 
GU046378 
Monte Grosso, Elba Island, Italy 42°51’N,10°24’E 2 GRO DQ3050765 DQ3050615 


terior probabilities of 2 99%. In light of these 
results, we present only analyses based on the 
two genes combined. We used PAUP* 4.0810 
to calculate base frequencies and to test for 
base frequency homogeneity across taxa (x2 
test). To detect saturation of substitutions, we 
plotted the absolute number of transversions 
and transitions against uncorrected-p genetic 
distances. This was done for each gene sepa- 
rately and for the two genes combined. For the 
COI dataset, saturation plots were derived for 
third codon positions only. 

Data was analyzed phylogenetically with S. 
similis (GenBank Accession Nos. DQ305085 
and DQ305059 for COI and 168, respectively) 
as outgroup. This species is sister to S. guidoni 
in the genus phylogeny (Ketmaier et al., 2006). 
We used MODELTEST (Posada & Crandall, 
1998) to determine the model of sequence 
evolution that best fit our data. We then used 
the MODELTEST output to calculate Maximum 
Likelinood (ML) genetic distances for Neighbor- 
Joining (NJ) analyses, as implemented in PAUP. 
We employed the same model of sequence evo- 
lution for Bayesian searches. These were run 


in MrBAYES 3.0b4 (Ronquist & Huelsenbeck, 
2003) allowing site-specific rate variation parti- 
tioned by gene and, for COI, by codon position. 
MrBAYES was run for two million generations 
with a sampling frequency of 100 generations. 
We ran one cold and three heated Markov 
chains and two independent runs. To assess 
if the Markov chains had reached stationarity, 
we plotted the likelihood scores of the sampled 
trees against generation time. Trees generated 
before stationarity were discarded as burn-in 
(first 10% of the sampled trees) and posterior 
probability values for each node were calculated 
on the basis of the remaining 90% of sampled 
trees. These trees were then used to construct a 
50% majority rule consensus tree in PAUP. The 
robustness of the NJ tree was tested by 1,000 
bootstrap replicates in PAUP. Phylogenetic tree 
topologies generated with different phylogenetic 
methods and competing phylogenetic hypoth- 
eses were evaluated statistically with the unbi- 
ased tree selection test (AU) as implemented 
in CONSEL (Shimodaira & Hasegawa, 2001). 
We compared tree topologies simultaneously 
(Shimodaira & Hasegawa, 1999). 
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The molecular clock hypothesis was tested 
with the likelihood ratio test (LRT; Goldman, 
1993), which compares the log-likelihood of 
ML trees with and without assuming a mo- 
lecular clock. Time estimates were calculated 
according to the formula T= Dy, /2r, where Dy 
is the maximum likelihood genetic distance 
and ris the mutation rate. Dy, distances were 
calculated for the COI and 16S partitions 
separately using the substitution model pa- 
rameters selected for each data partition with 
MODELTEST (Posada & Crandall, 1998). We 
used previous published rates for mollusks to 
derive divergence times; the same rates were 
adopted in Ketmaier et al. (2006). Thus, time 
estimates between that study and the current 
one are perfectly comparable. Rates vary be- 
tween 0.67 and 1.21% per million years for COI 
and between 0.5 and 0.6% per million years 
for mitochondrial ribosomal genes (Ozawa & 
Okamoto, 1993; Rumbak et al., 1994; Marko, 
2002). 

We used the Nested Clade Phylogeographi- 
cal Analysis (NCPA) (Templeton, 1998) to 
test for association between genealogy and 
geography in the dataset and to infer popula- 
tion processes. This method has been recently 
criticized, because it seems to be prone to a 
high rate of false positives (e.g. Knowles, 2008; 
Petit, 2008). However, NCPA has still a number 
of advantages compared to alternative meth- 
ods (Templeton, 2009). In addition, we based 
our main conclusions on multiple independent 
analyses (i.e. phylogenetic and population 
genetics methods) to complement the NCPA 
inferences. We used the software ANeCA 
(Panchal & Beaumont, 2007), which automates 
the complex NCPA and renders manual proce- 
dures (i.e. the nesting design) in the analysis 
unnecessary. The package includes TCS, 
which implements the statistical parsimony 
procedure (Clement et al., 2000), and GEO- 
DIS (Posada et al., 2000) to respectively 
derive the haplotype network and to calculate 
distances among haplotypes and clades and 
their statistical significance. We included four 
additional locations in the GEODIS analysis 
(Monte Castello, Elba Island, 42°47’N, 10°22’E: 
Col de Teghime, Corsica, 42°40’N, 9°24’E; 
Sabara, Corsica, 42°24’N, 9°10°E; Punta Giglio, 
Sardinia, 40°34’N, 8°12’E; indicated by black 
dots in Fig. 1). These are localities where the 
species had been found in the past but not dur- 
ing our sampling; including this information is 
important to detect any sampling inadequacies 
(Panchal & Beaumont, 2007). The statistical 


distribution of distances was determined by 
re-calculating distances after 10,000 random 
permutations (i.e., clades against sampling 
localities). ANeCA automatically runs a slightly 
different version of the inference keys, dated 
14 July 2004. 

We used ARLEQUIN 3.0 (Excoffier et al., 
2005) to calculate the following parameters of 
genetic diversity for each population: haplotype 
diversity (h), mean number of pairwise differ- 
ences between all pairs of haplotypes (77) and 
the nucleotide diversity (77,,). Levels of genetic 
diversity were tested by a hierarchical analysis 
of molecular variance (AMOVA) (Excoffier et 
al., 1992) in ARLEQUIN with 1,000 replications. 
For this analysis we grouped populations (1) by 
island of provenance, and (2) by island but with 
Corsican populations split into two groups as 
suggested by the phylogeographic and NCPA 
results. We also used ARLEQUIN to calculate 
pairwise Fsy values for all pairs of populations; 
statistical significance of Fst values was as- 
sessed by 1,000 permutations with sequential 
Bonferroni correction for multiple tests. 


RESULTS 
Sequence Variation 


We sequenced 1,143 bp of mtDNA (675 
bp for the COI gene and 468 bp for the 16S 
rRNA gene) for each individual included in the 
study. We observed only few gaps in the final 
alignment all in the 16S gene and confined to 
outgroup vs. ingroup comparisons. Levels of 
sequence variability were slightly higher in COI 
than in 16S; there were 50 (7.40%) polymorphic 
sites in COI and 33 (7.05%) in 16S; parsi- 
mony informative sites were 35 (5.18%) and 10 
(2.13%) for the two genes. As expected, most of 
the variation resided in COI 3rd codon positions; 
43 (19.11%) sites were variable, and of these 
31 (13.7%) were also parsimony informative. 
The combined data set had 83 (7.26%) and 45 
(3.93%) variable and parsimony informative 
sites. Sequences were generally A+T rich and 
anti-G biased, a typical pattern for mitochondrial 
genomes. A+T percentages ranged from 57.5% 
in COI 1st codon positions to 88.9% in COI 3rd 
codon positions; G percentages varied from 
5.3% in COI 3rd codon positions to 27.9% in COI 
1st codon positions. A x2 test of homogeneity 
of base frequencies among taxa did not reveal 
any significant differences, irrespective of the 
data partition tested (both genes combined, 
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TABLE 2. Haplotype frequencies of all individuals of Solatopupa guidoni sequenced for the study. For 
population codes see Fig. 1. Numbers in brackets are sample sizes. 


Sardinia 
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each gene separately and COI 1st, 2nd, and 3rd 
codon positions separately). Visual inspection 
of saturation plots (not shown) revealed that 
saturation is not present in the data set. 


Haplotype Frequencies and Genetic Variability 
of Populations 


The 60 sequences yielded 21 haplotypes 
(Table 2). These were population-specific (i.e., 
not shared among populations) and differed 
by 1 to 23 substitutions. COR, CAP, OLE, FLO 
and SIS had more than one haplotype with 
a maximum number of five in CAP (Table 2). 
CCC, SBR, FRA and GRO were fixed for single 
haplotypes. In the other populations, haplotype 
diversity (h) ranged between 0.222 + 0.166 
(OLE) and 1.000 + 0.126 (CAP); the mean 
number of pairwise differences between all pairs 
of haplotypes (1) and nucleotide diversity (t1,) 


Populations 


Corsica Elba Is. 


FRA(3) CAP (5) OLE (9) FLO (10) SIS(8) GRO (2) 


—_— j= =k = od 


- N A = 


was lowest in COR (0.375 + 0.376 and 0.0003 
+ 0.0003) and highest in OLE (4.666 + 2.523 
and 0.004 + 0.002). 


Phylogenetic Analyses and Divergence Times 


Figure 2 shows the haplotype Bayesian tree 
obtained for the combined data set using the 
GTR + | model (proportion of invariable sites 
= 0.876, equal distribution of rates at variable 
sites) and summarizes the results of the NJ 
searches. Both trees were statistically indis- 
tinguishable with the AU test (P = 0.641) and 
identified two major haplotype groups. A first 
group received good bootstrap support (89%) 
but fair posterior probability (89%); this group 
included haplotypes from Sardinia, Corsica, 
and Elba Island. Relationships within this clade 
were generally supported. Haplotypes from two 
Corsican localities (FLO and OLE) did not form 
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FIG. 2. Bayesian haplotype phylogram based on the GTR+I model of sequence evolution. 
Only statistical supports > 75% for both the Bayesian (2,000,000 generations) and NJ 
searches (1000 bootstrap replicates) are reported on branches (first and second values). 
Population codes are as in Table 1 while haplotype numbers are as in Table 2. A, B, C, D and 
E identify nodes for which estimates of divergence times have been calculated (Table 3). 


monophyletic clusters. Rather, haplotype 10 
(FLO) grouped with one of the two haplotypes 
found at OLE (16) plus the haplotype sampled 
on Elba Island (GRO). The second haplotype 
from OLE (15) was sister to haplotype 6 found 


on Sardinia. The second group only included 
Corsican haplotypes. Relationships within this 
clade were poorly resolved, with only a low 
support (81%) for haplotypes from CAP as 
monophyletic in the Bayesian analysis. 
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TABLE 3. Inferred divergence times (Myr) based on the COI and 16S data partitions and on a variety 
of published rates for mollusks for supported splits in the Solatopupa guidoni haplotype tree of Fig. 2 
(nodes A-E}. Maximum Likelihood distances (DuL) are based on the models of sequence evolution that 
best fitted the two data partitions (models selected with MODELTEST). *Rates from Ozawa & Okamoto, 
1993; Rumbak et al., 1994; and Marko, 2002. “Substitution models are GTR + | and F81 for COI and 
16S data partitions, respectively. 


Data partitions/ COI 16S 
Rates* (%/Myr) 0.67—1.08 0.5-0.6 
Node Dy Myr DuL” Myr 
A 0.007 + 0.002 0.52 +0.14—-0.32 + 0.009 0.0005 + 0.0005 0.05 + 0.05-0.04 + 0.04 
B 0.01+0.001 0.75+0.07-0.46+0.046 0.0012+0.001 0.12+0.100-0.1 + 0.083 
C 0.001 + 0.0006 0.07 + 0.044—0.05 + 0.027 0.000 - 
D 0.015 1.12-0.69 0.009 0.9-0.75 
E 0.034 + 0.001 2.53 + 0.07-1.57 + 0.046 0.004 + 0.0004 0.4 + 0.03-0.33 + 0.03 
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FIG. 3. Haplotype network and Nested Clade Phylogeographical Analysis of Solatopupa guidoni 
populations derived from COI + 16S data. The relative size of the circles is proportional to the number 
of individuals carrying that particular haplotype; different shadings indicate the geographic origin of 
haplotypes. Coding of populations and numbering of haplotypes are as in Tables 1 and 2. Numbers in 
parentheses indicate how many individuals carried that particular haplotype. Black dots are missing 
haplotype. Haplotypes are always one mutational step away from each other if not indicated otherwise 
(i.e. when the number of steps was higher than 4 numbers in italics are the mutational steps necessary 
to connect haplotypes; missing haplotypes are not represented in these particular circumstances). The 
Figure also shows the nesting design as obtained with ANeCA. Clades are designed as x-y, where x 
represents the nesting level (from the haplotype to the 4th level but haplotype and 1-step levels are not 
shown) and y is the number assigned to that particular clade. 
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TABLE 4. Nested clade analysis of the geographical distribution of Solatopupa guidoni haplotypes. Only 
the six clades with significant GEODIS results are shown (see P value of the x2 test; third column). 
Clade names and nesting design are as in Fig. 3; the geographical distribution of clades is also shown. 
The last two columns show the sequence followed in the inference key implemented in ANeCA and 
the resulting biological inference. 


Clade Distribution ad Inference key sequence Biological inference 
2-11 FRA, CAP (Corsica) 0.050 1,2, tip-interior status cannot be Inconclusive outcome 
determined 
3-1 SBR, FRA, CAP (Corsica) 0.005 1,19, no Allopatric fragmenta- 
tion 
3-3 ‘OLE, FLO, SIS, GRO (Cor- 0:000 "172511, 12ane Contiguous range 
sica, Elba Is.) expansion 
4-1 COR, SBR, FRA, CAP (Cor- 0.000 1,2, tip-interior status cannot be Inconclusive outcome 
sica) determined 
4-2 CCC, OLE, FLO, SIS, GRO 0.000 1,2,3,4,9, no Allopatric fragmenta- 
(Sardinia, Corsica, Elba Is.) tion 
Total All populations (Sardinia, 0.000 1,2, tip-interior status cannot be Inconclusive outcome 
Corsica, Elba Is.) determined 


Because the Corsican populations did not 
form a monophyletic group, we tested this al- 
ternative hypothesis against the unconstrained 
topology presented in Figure 2. Yet, the AU 
test rejected the monophyly of the Corsican 
populations (P = 0.008). 

The LRT test conducted on trees with and 
without the molecular clock alternatively relaxed 
and enforced could not reject the hypothesis of 
homogeneous rates along the branches of our 
phylogeny [-InLnoclock = 2114.27; -InLclock = 


freedom = 31.41]. We could therefore calculate 
divergence times for labeled nodes in the tree 
of Figure 2 (Ato E); the results of these calcula- 
tions are presented in Table 3. These estimates 
have to be considered cautiously, as not all the 
above nodes were unambiguously supported by 
our phylogenetic searches though the network 
analysis retrieved almost identical groupings. 
Time estimates based on COI ranged in age 
from 0.05 + 0.027 Myr (node C) to 2.53 + 0.07 
Myr (node E); those based on 16S varied be- 


tween 0.04 + 0.04 Myr (node A) and 0.4 + 0.03 
Myr (node E). 


2115.13; 2AL = 1.72; degrees of freedom = 20; 
critical x2 value (P = 0.05) for 20 degrees of 


TABLE 5. Results from the analysis of molecular variance (AMOVA). Populations have been grouped 
alternatively into three (scheme !) or four groups (scheme II). Scheme | corresponds to the three is- 
lands (Sardinia, Corsica and Elba Island) sampled for the study while in scheme II we split Corsican 
haplotypes into two groups (COR, SBR, FRA, CAP and OLE, FLO, SIS, respectively), as suggested 
by phylogenetic and network analyses (Figs. 2, 3). 


Analysis Sum of squares Variance 
structure Source of variation (d.f.) component P Fixation index % of variation 
scheme | Among groups 2.68 (2) Va = 0.062 0.040 Ọecr= 0.113 11.35 
Among populations 15101-(6) Vb = 0.310 0.000 Ọsc = 0.636 56.44 
Within populations 9.02 (31) Vc=0.180 0.000 Do; = 0.677 32.21 
scheme II Among groups 6:12 (3) Va = 0.023. 0.000 Ọcr= 0.047 8.94 
Among populations 11.56 (5) Vb = 0.305 0.000 Ọsc= 0.633 63.32 
Within populations 9.02 (51) Vc=0.176 0.000 cz = 0.650 27.74 
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TABLE 6. Fg; values among the Solatopupa guidoni populations included in the study. Bold 
numbers indicate statistical significance after sequential Bonferroni correction for multiple 
comparisons (significance assessed by 1,000 permutations in ARLEQUIN 3.01). 


Sardinia Corsica Elba Is. 
Population 1 2 3 4 5 6 T 8 9 
1-CCC - 
2-COR 0.798 Ta 
3-SBR 1.000 0.809 - 
4-FRA 1.000 0.798 1.000 - 
5-CAP 0.387 0.519 0.452 0.387 - 
6-OLE 0.828 0.765 0.843 0.828 0.469 - 
7-FLO 0.711 0.703 0.733 0.711 0.377 0.696 - 
8-SIS 0.479 0.563 0.522 0.479 0.138 0.524 0.449 = 
9-GRO 1.000 0.784 1.000 1.000 0.285 0.809 0.681 0.415 - 


Nested Clade Phylogeographical Analysis, 
AMOVA and Fes7 values 


Figure 3 illustrates the network obtained with 
TCS (Clement et al., 2000) and shows the 
nesting design with four nested levels. The 


Local diversification 
including colonization of Sardinia 
and Elba island 
(0.04-2.5 Myr) 


pa Northern coastal Corsican populations 


Ei Central inland Corsican populations 


A 
37 Ea 


pattern of relationships in the network was 
generally congruent with that depicted by the 
phylogenetic tree of Figure 2. Clades 4-1 and 
4-2 corresponded to the two major groups of 
haplotyes found in the phylogenetic searches. 
These clades were separated by a deep phy- 


GRO 
Split between the S. guidoni lineage 
and its continental ancestor 
a. My) 


Split of S. guidoni 
into coastal and mainiand lineages 
(0.9-4.1 Myr) 


FIG. 4. Schematic of proposed phylogeographic history of Solatopupa guidoni with patterns and timing 
of the major events of differentiation within the species as inferred from phylogenetic and nested clade 
analyses combined with molecular clock estimates. See text for detailed discussion. 
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logeographic break of 20 substitutions. The 
two haplotypes from OLE (15 and 16) were 
not linked to each other. Overall, the network 
analysis confirmed that the Corsican haplo- 
types do not form a monophyletic entity and can 
be split into two major sub-groups. The nesting 
procedure produced eight clades that could be 
tested for geographical association, and six 
of them resulted in significant permutational 
contingency tests indicating a non-random geo- 
graphical distribution of haplotypes (Table 4). 
The inferences derived from Templeton’s keys 
suggested allopatric fragmentation for clades 
3-1 and 4-2 and contiguous range expansion 
for clade 3-3. In three circumstances (2-11, 
4-1, and for the total cladogram) the inference 
keys gave inconclusive outcomes because of 
the unresolved status of tip-interior position of 
the clades. 

Hierarchical analysis of molecular variance 
(AMOVA) is reported in Table 5. Pooling Cor- 
sican populations in a single group or splitting 
them into two groups did not change the main 
outcomes of the analysis. AMOVA showed that 
most genetic heterogeneity was apportioned 
among populations within groups and within 
populations. Together these two levels explain 
between 87% and 92% of the overall variation, 
depending on the hierarchical scheme adopted 
(lor ll; Table 5). These results, along with those 
of the pairwise Fs analysis reported in Table 6, 
indicate a general high level of genetic differen- 
tiation. The vast majority of pairwise Fsy values 
were over 40% and statistically significant after 
sequential Bonferroni correction. 


DISCUSSION 


The main findings of this study can be sum- 
marized as follows: S. guidoni is a genetically 
structured species that can be broadly divided 
into two lineages. One of these lineages is 
limited to central inland Corsica, while the 
second includes northern coastal Corsica sites 
plus Sardinia and Elba populations. Patterns of 
relationships suggest that differentiation within 
Corsica predated between-island divergence. 
This translates into a lack of reciprocal mono- 
phyly of lineages occurring on the three land- 
masses. While Sardinia and Elba both harbor 
monophyletic (monotypic) lineages (but a very 
limited number of individuals were analyzed for 
these locations), they are nested in one of the 
Corsican clades. 


Levels of Divergence and Population Structure 


Numerous studies have revealed a broad 
range of intraspecific mtDNA differentiation 
in terrestrial mollusks (2.9-30% of sequence 
divergence) (Thomaz et al., 1996; Pfenninger & 
Magnin, 2001; Davison, 2002; Holland & Had- 
field, 2002; Uit de Weerd et al., 2005). We found 
an average level of sequence divergence of 3.6 
+ 0.2% between the two alleged haplo-groups 
of S. guidoni. This is close to the lower bound 
of the aforementioned range of mtDNA varia- 
tion. These results suggest that processes of 
differentiation within the species are relatively 
recent. Our molecular clock estimates below 
support this hypothesis. In addition, Giusti 
(1977) stressed that there are no clear differ- 
ences in shell morphology or internal anatomy 
that could help distinguish Sardinian, Corsican, 
and Elba Island populations from one another. 
Thus, neither genetics nor morphology (Giusti, 
1977) support the need for taxonomic revision 
of the species, contrary to the claim of Boato 
et al. (1985). 

On the other hand, our data indicate a gen- 
eral lack of gene flow among populations of S. 
guidoni, as attested by haplotype frequencies, 
AMOVA results and generally high and signifi- 
cant pairwise Fs; values. AMOVA showed that 
populations were significantly differentiated 
at all three levels of the hierarchical design, 
with a large fraction of variation apportioned 
among populations, whatever hierarchical 
scheme was adopted in the analysis (Table 
5). These results are in line with allozyme data 
that demonstrated a substantial lack of gene 
flow in the species, with populations bearing 
a considerable number of private and fixed 
alternative alleles and with an average number 
of migrants per generation among locations of 
0.036 (Boato, 1988). Concordant patterns in 
markers with different characteristics (nuclear 
vs. mitochondrial) suggest the same underlying 
evolutionary scenario. These results also fit 
well with some general characteristics of land 
snails. Populations of these organisms are 
often divided into small breeding units (demes) 
that are partially, if not completely, isolated from 
one another by habitat discontinuities. Move- 
ments of individuals between demes are very 
rare, since distance between nearby demes 
is usually larger than the average dispersal 
range of individuals (Fiorentino et al., 2009, 
and references therein). It is therefore not 
surprising that most of the molecular papers 


92 KETMAIER ET AL. 


published so far on terrestrial mollusks have 
revealed considerable population structur- 
ing and geographic variation (reviewed in 
Davison, 2002). Stochastic processes (i.e., 
genetic drift and founder effect) have been 
identified as the prominent factors shaping 
genetic traits of populations. Coalescence 
of haplotypes is another process that might 
account for population division in terrestrial 
mollusks. Their populations are typically ar- 
ranged in a stepping-stone fashion and divided 
into discrete demes of large size (hundreds or 
thousands of individuals). Their high density 
and the limited migration of individuals be- 
tween them theoretically translate into much 
longer persistence (time to common ancestor) 
of mitochondrial haplotypes with retention of 
ancient lineages and thus into an increase in 
overall divergence (Thomaz et al., 1996). It 
must be noted, however, that population size 
in Solatopupa tends to be small compared to 
that of other snails and occurrence of demes 
is scattered (Boato, 1988). These seem to be 
common traits in rock-dwelling species bound 
to calcareous substrata (Fiorentino et al., 
2009). Moreover, rocks of different geological 
nature surround the limestone occupied by 
the species constituting a major impediment 
to dispersal. Assuming, contrary to what is 
commonly reported for other land snails, that S. 
guidoni populations are arranged in scattered 
demes of small size, we argue that time to 
coalescence may be relatively short, explaining 
the rather shallow level of divergence. 


Historical Biogeography 


A number of molecular studies on a taxo- 
nomically diverse array of organisms endemic 
to the peri-Tyrrhenian area have been carried 
out in the last decade. The organisms include 
newts, subterranean terrestrial and aquatic 
isopods, stoneflies, and flightless beetles 
(Cobolli Sbordoni et al., 1995; Caccone et al., 
1997; Ketmaier et al., 1999, 2000, 2003a, b; 
Fochetti et al., 2004). All these organisms are 
found in Sardinia and Corsica and some of 
them also occur in the Tuscan archipelago or 
coastal areas of Tuscany, that is, the terrestrial 
and aquatic isopod genera Oritoniscus and 
Stenasellus (Cobolli Sbordoni et al., 1995; 
Ketmaier et al., 1999, 2000, 2003a). In all 
circumstances, phylogenetic analyses were 
concordant in sorting Sardinian and Corsican 
taxa into distinct groups separated by large 
mutational distances. Tuscan representatives, 


if any, were always more closely related to the 
Corsican than to the Sardinian taxa (Cobolli 
Sbordoni et al., 1995; Ketmaier et al., 1999, 
2000, 2003a). Interestingly, a major phylogeo- 
graphic break between Corsica and Sardinia 
and a closer affinity between populations from 
Corsica and the Tuscan archipelago has also 
been found in plants belonging to the genus 
Borago (Selvi et al., 2006). Conversely, we did 
not find reciprocally monophyletic lineages in 
S. guidoni on the three islands. 

Given the lack of fossil evidence for S. gui- 
doni, which prevents the use of tree-based 
estimates of time divergences (Sanderson, 
1997, 2002, 2003), we could only rely on rates 
published for mollusks. It is important to stress 
that mtDNA substitution rates are a notoriously 
difficult issue in gastropods because of their 
extensively documented accelerated rates of 
evolution (Chiba, 1999; Pinceel et al., 2005). 
However, we tend to reject this hypothesis 
for S. guidoni on the grounds of the avail- 
able allozyme data (Boato, 1988). MtDNA 
and allozymes support concordant patterns 
of relationships as shown here. According to 
Pinceel et al. (2005), this is good evidence that 
mt- and nuclear DNA rates are not strongly de- 
coupled with mtDNA evolving much faster than 
nucDNA. If they were, we would have observed 
divergence in the mtDNA but sharing of alleles 
at the allozyme loci across locations, which is 
not the case here as different populations are 
characterized by a high percentage of private 
alleles (12.5%; Boato, 1988). Additionally, time 
estimates for the statistically supported nodes 
in our phylogeny are all relatively recent, not 
older than 2.5 Myr. While these are only crude 
estimates of minimum time since divergence 
and are based on a not fully supported phylog- 
eny, they are in line with the time frame of genus 
cladogenesis, which took place between the 
middle Miocene and the middle-late Pliocene 
(Ketmaier et al., 2006). 

Previous allozyme and DNA sequence data 
had already suggested that Corsican popula- 
tions might be divided into a central inland and 
a northern coastal clade (Boato, 1988; Ketmaier 
et al., 2006). These would correspond to clades 
4-1 and 4-2 of our study (Fig. 3). Sardinian and 
Elba Island haplotypes are linked to the north- 
ern coastal clade. Boato (1988) found a sister 
relationship between a population sampled 
close to OLE (Col de Teghine, black dot in Fig. 
1) and CCC (same locality as in the present 
study), but no samples from Elba Island were 
included in that study. In our previous paper, 
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we retrieved a sister relationship between CCC 
and GRO, but only analyzed two mainland Cor- 
sican populations (FRA and SBR) (Ketmaier et 
al., 2006). The occurrence of the species on 
Sardinia and Elba Island can be explained in 
the light of the paleogeographic history of these 
landmasses and their connections to Corsica. 
According to Steininger & Rdg! (1984), Corsica 
and Sardinia became permanently isolated in 
the Pliocene (5 to 3 Myr BP). More recent data, 
however, suggests that this is an over-simplified 
scenario and that multiple relatively short-lived 
connections between the two islands existed 
during the Quaternary, up to 20,000 years ago 
at each glacial peak, as a consequence of the 
worldwide drop in sea level (Meulenkamp & 
Sissingh, 2003; Pascucci, per.com.). Our time 
estimates for the split between the Sardinian 
population and its closest Corsican relative 
(haplotype 15, OLE) range from 0.7 to 1.12 
Myr BP, depending on the gene used and the 
substitution rate adopted (Table 3). Although 
these datings span a relatively broad time 
scale, they overlap remarkably with the GUnz 
glaciation (West, 1977). They suggest that a 
(possibly single) dispersal event of the species 
from Corsica to Sardinia took place at that time, 
but not during more recent Quaternary connec- 
tions between the two islands. This hypothesis 
is further supported by the fact that in Sardinia 
the species is limited to two nearby locations 
(Fig. 1), despite the extensive occurrence of 
limestone throughout the island. The haplotype 
sampled on Elba Island is deeply embedded 
within Corsican haplotypes in the phylogenetic 
tree (Fig. 2), while it has a tip position in the 
haplotype network (Fig. 3). According to our 
time estimates, differentiation within this haplo- 
group occurred in the Pleistocene, between 
0.04 and 0.52 Myr BP. Anti-clockwise rotation 
of the Sardinia-Corsica microplate caused up- 
lifting of coastal Tuscany, including the Tuscan 
archipelago, which has been repeatedly con- 
nected to Corsica as well as to coastal Tuscany 
during the Quaternary (Lipparini, 1976; Cherchi 
& Montadert, 1982; Burgassi et al., 1983; 
Pascucci et al., 1999). As a consequence, the 
archipelago hosts a fauna transitional between 
the two areas. The incidence of continental 
Tuscan and Corsican species in the Tuscan 
archipelago is directly related to the present- 
day distance of the different islands from the 
two source areas, with Elba Island playing the 
role of stepping stone (Dapporto et al., 2007, 
and references therein). A close phyletic rela- 
tionship between Corsican and insular Tuscan 


species has frequently been found for species 
with poor dispersal capacities, such as stygobi- 
ont taxa (Ketmaier et al., 2003a). This led to the 
recognition of a Corsica-Sardinian stygofaunal 
province that also includes the insular areas 
of Tuscany (Pesce, 1985). We conclude that 
the most parsimonious scenario to account 
for the exclusive occurrence of the species on 
Elba Island (and not on any other island of the 
Tuscan Archipelago) is based on dispersal. 
Our time estimate for this event suggests that 
the species reached Elba Island more recently 
than Sardinia. 


The Phylogeographic Scenario 


The nested clade analysis suggests a dis- 
persal event between Corsica and Elba Island 
with contiguous range expansion (clade 3-3; 
Fig. 3). This implies movements of individuals 
as connections between landmasses became 
available. For clades 3-1 (central inland Corsi- 
can lineage) and 4-2 (Sardinia, coastal northern 
Corsica and Elba Island), the analysis indicated 
allopatric fragmentation as the most likely bio- 
logical inference. A combination of vicariance 
and dispersal events would therefore most ef- 
fectively explain the current distribution of the 
species. Phylogeographic relationships within 
S. guidoni are evidently complex; Figure 4 
summarizes the following discussion. Based on 
our phylogenetic, nested clade and molecular 
clock analyses we propose that immediately 
after the split from its continental sister taxon, 
S. guidoni underwent further division into a 
coastal northern and central inland Corsican 
lineage. Diversification within the latter group 
is presumably allopatric and associated with 
the scattered distribution of suitable calcareous 
rocks among geologically different substrata. 
Relationships within the coastal lineage are 
less clear. A first split separated FLO from 
OLE (haplotype 15, the most abundant in the 
latter population). The occurrence of the spe- 
cies In Sardinia is difficult to explain because 
of the large geographical gap separating the 
Sardinian population from its closest extant 
population (OLE). In addition, S. guidoni has 
never been found in southern Corsica, despite 
extensive surveys in the only area with lime- 
stone (Bonifacio; shaded in grey in Fig. 1). We 
cannot exclude that S. guidoni had a broader 
range, occurred there sometime in the past 
and then became extinct; unfortunately, this is 
a hypothesis that currently cannot be tested. 
According to the haplotype network (Fig. 3), 


94 KETMAIER ET AL. 


FLO gave rise to SIS, probably through a single 
founder event, as supported by the monophyly 
of SIS (Figs. 2, 3). Nearly coeval with this was 
the origin of the lineage that includes the sec- 
ond haplotype found at OLE (haplotype 16) and 
the haplotype found on Elba Island. The occur- 
rence of divergent haplotypes at OLE implies 
admixture between presumably allopatrically 
evolved haplotypes (Avise, 2000). The only 
known population in the area that we were not 
able to sample is Col de Teghine (Fig. 1). This 
population (but not OLE) was included in the 
allozymic survey of Boato (1988) and it showed 
affinity with FLO, while bearing private alleles. 
We tentatively conclude that at OLE a more 
divergent haplotype (15) met haplotype 16 of 
more recent origin. 

This study reveals that within-island di- 
vergence in S. guidoni is slightly older than 
between-island differentiation and that Corsica 
(the most intensively sampled island) does not 
harbor a monophyletic lineage. While very little 
can be derived for Sardinia and Elba Island be- 
cause of the small sample sizes, for Corsica our 
results do not support a one island-one clade 
association, even in an organism with theoreti- 
cal poor dispersal capacity and distributed over 
biogeographically well-known areas. 
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